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Abstract 

We use a modification of the parameterization method to study invariant manifolds for difference 
equations. We establish existence, regularity, smooth dependence on parameters and study several 
singular limits, even if the difference equations do not define a dynamical system. This method also 
leads to efficient algorithms that we present with their implementations. The manifolds we consider 
include not only the classical strong stable and unstable manifolds but also manifolds associated to 
non-resonant spaces. 

When the difference equations are the Euler-Lagrange equations of a discrete variational we 
present sharper results. Note that, if the Legendre condition fails, the Euler-Lagrange equations can 
not be treated as a dynamical system. If the Legendre condition becomes singular, the dynamical 
system may be singular while the difference equation remains regular. We present numerical appli- 
cations to several examples in the physics literature: the Frenkcl-Kontorova model with long-range 
interactions and the Heisenberg model of spin chains with a perturbation. We also present extensions 
to finite diffcrcntiablc difference equations. 

1 Introduction 

1.1 Difference equations 

In this paper we generalize some notions that have played an important role in dynamics, namely 
invariant manifolds, to the more general context of implicit difference equations. We also present 
algorithms to compute this object and apply them to several problems in the physics literature. 
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Consider a smooth manifold M of dimension d. Let Z : M N+1 — > M d be an analytic function. The 
function Z represents the following difference equation of order N 



we can consider Q as defining 6k+N m terms of the other variables, so that one can transform ([T]) into a 
recurrence of the form Q. However, there are instances in which it is impossible to do so and therefore 
problem is more general than Q. Recurrences have been studied mainly as dynamical systems and 
they have a rich geometric study. 

The goal of this paper is to show that some familiar constructions in the theory of dynamical systems 
have rather satisfactory generalizations in the context of implicit difference equations. We will use the 
so-called parameterization method [CFdlL03a, CFdlL03b, CFdlLQ5] to show existence, regularity and 
smooth dependence on parameters of these invariant manifolds. Note that other methods such as the 
graph transform, which depend very much on the dynamical formulation do not seem to generalize to 
the context of implicit difference equations. 

The Legendre transform -if it exists- makes a difference equation into a dynamical system. In 
some situations, a family of Legendre transforms exist but becomes singular. We study some of these 
singular limits, in which some invariant manifolds continue to exist through the singularity of this 
Legendre transformation and the implicit equation remains smooth, even if the dynamical system does 
not. 

We will also describe and implement some rather efficient numerical algorithms to compute these 
objects with high accuracy. We point out that some of the algorithms are novel, even in the dynamical 
system case since we study not only the classical stable (and strong stable) manifolds but also the weak 
stable manifolds as well as singular limits. We believe that this is interesting because it shows a new 
way of approximating stable and unstable manifolds. In order to show that the method is robust we 
perform some explicit calculations that can be found in the code that supplements this paper. 

Remark 1.1. We note that the relation between difference equations in implicit form ([I]) and the 
explicit recurrence is similar to the relation of Differential-Algebraic Equations (DAE) of the form 
Z (y,y r , ■ ■ ■ ,y^ n ') = and explicit Differential equations j/ n ' = F (y,... , It seems that the 

methods presented here can be extended to DAE, but one needs extra technicalities. We hope to come 
back to this question. For more information on DAE, see [KM06j. 




(1) 



The solutions of the difference equation (|T 
k > 0. A recurrence 



are sequences (Ok)k>o which satisfy ([!]), for all values of 
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1.2 Some examples and applications 

1.2.1 Variational principles 

One important source of the problems of the form ([!]), to which we pay special attention, is discrete vari- 
ational problems (which appear in physics, economics, dynamic programming, etc.). A mathematical 
review of these discrete variational problems can be found in [Ves91, GolOl]. 

Let S : M N+l — > M be analytic. When studying variational problems one is interested -among 
other things- in solutions to Euler-Lagrange equations Z (9 k , . . . , 9k+2N) = 0, with 

N 

Z(8 , . . . , 8 2 n) = ^ djS (8N-j, &N-j+i, ■ ■ ■ , ^2N-j) = 0, (3) 

3=0 

which are of the form ([!]) and of order 2N. The Euler-Lagrange equations appear when one finds 
critical points of the formal variational principle based on J {9) = Ylk & (@k> @k+i, ■ ■ ■ , 9k+N)- Indeed, 
^ is formally just ^ J(O) = 0, for all k. 

As we mentioned, there are well known conditions (Legendre conditions) which allow to transform 
the Euler-Lagrange equations into recurrences, often called discrete Hamilton equations or twist maps. 
In situations where Legendre conditions fail -or become singular- the Euler-Lagrange equations cannot 
be transformed into Hamiltonian equations. See the examples below and those treated in more detail 
in Section [HI 

Even if we will not explore it in detail, we note that the invariant manifolds considered here have 
applications to dynamic programming. We plan to come back to these issues in future investigations. 

1.2.2 The Heisenberg XY model of magnetism 

There are many examples of Lagrangian systems that can not be transformed in dynamical systems. 
Consider, for instance, the Heisenberg XY model. In this model, one is interested in solutions of 

sin(0 fc+ i - 9 k ) + sin(6» fc _i - 9 k ) - e sin 9 k = 0, (4) 

where each k is an angle that represents the spin state of a particle in position k £ 7L. The parameter 
e corresponds to the strength of an external magnetic field. The equation Q appears as the Euler- 
Lagrange equation of an energy functional J{9) = Sfc cos (^fc+i — ®k) + ^cos{9 k ) and has Lagrangian 
function S{9q,9\) = cos(6*i — 9q) + ecos(6>o)- 

The dynamical interpretation of Q is problematic because in order to get 9 k+ \ in terms of 9 k and 
0fc-ij w e need to have 

\sm(9 k - 9k-i) + esm9 k \ < 1. 

This condition is not invariant under the dynamics and having it for one value of k does not guarantee 
to have it for others. However, 9 k = is a solution of Q. In fact, there could exist many solutions, 
defined for k > 0, that converge to the fixed point. So, it is interesting and useful to identify these 
solutions and understand their geometry. We will present algorithms illustrating our general results in 
this model in Section [6. 4[ 
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1.2.3 The Frenkel-Kontorova model with non-nearest interaction 

The Frenkel-Kontorova model was introduced to describe dislocations |FK39] . but it has also found 
interest in the description of deposition |AL83l |BK04| . One simplified version of the model (more general 
versions will be discussed in Section 6.3) leads to the study of solutions of non-nearest interactions. For 
instance, we have 

e (9 k+2 + G k -2 - 26 k ) + (9 k+1 + 9 k -x - 26 k ) + V'(6 k ) = 0. (5) 

Note that for e = 0, equation ^ can be transformed into a second order recurrence; i.e. a dynamical 
system in 2 dimensions. Indeed, this 2-D map is the famous standard map, [Mat93"l iGolOlj . whereas for 
e 7^ - no matter how small - one is lead to a dynamical system in 4 dimensions. 

It turns out that some terms in this 4— dimensional system blow up as e — > 0. Hence, the perturba- 
tion introduced by the e term is singular in the Hamiltonian formalism. Nevertheless, we observe that 
the singularity appears only when we try to get 6 k+ 2 a s a function of k +i, k , 9 k -i, Q k -2- The equations 
([5]) themselves depend smoothly on parameters. 

Indeed, in section 5.1 we will show that the invariant manifolds of §5§ that we construct, are 
smooth across e = 0. We note that in the applications to solid state physics, the regime of small e is 
the physically relevant, one expects that there are interactions with even longer range which become 
smaller with the distance |CDFM07l fSuz71j . 



1.2.4 Dependence of parameters of invariant manifolds and Melnikov theory 

In many situations, the solutions of a difference equation change dramatically when parameters are 
introduced. In the case of dynamical systems the transverse intersection of the stable and unstable 
manifolds is associated with chaos, and gave rise to the famous horseshoe construction of Smale. The 
Poincare-Melnikov method is a widely used technique for detecting such intersections starting from a 
situation when the manifolds agree. One can assume that a system has pair of saddles and a degenerate 
heteroclinic or saddle connection between them. The classical Melnikov theory computes the rate at 
which the distance between the manifolds changes with a perturbation. In this context, it is important 
to understand dependence of parameters in invariant objects that appear in a Lagrangian setting. 
In particular, in the XY and XYZ models, dynamics is not longer useful and a purely variational 
formalism is needed. In this paper, we show that the variational theory is robust and there is smooth 
dependence on parameters, so our theory could lead to a variational formulation of Melnikov's theory. 
Previous work in this direction was done in |Tab95| ILom97] in the context of twist maps. 

1.3 Organization of the paper 

The method of analysis is inspired by the study in [CFdlL03a, CFdlL03b, CFdlL05] and is organized 
as follows. First, in Section [3] we will study the linearized problem and generalize the notion of charac- 
teristic polynomial, spectrum and spectral subspaces. They give necessary conditions so that one can 
even consider solving (ISJ) . 
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The Main Theorem 4.1 is stated in Subsection |4.1| After we make the choices of invariant spaces 
in the linear approximation, we will show that, provided that these spaces satisfy some non-resonance 
conditions -automatically satisfied by the (strong) stable spaces considered in the classical literature- 
the solutions of the linear problem lead to solutions of the full problem. 

The main tool is some appropriate implicit function theorems in Banach spaces. These implicit 
function theorems some times require that we consider approximate solutions of order higher than the 
first. In Appendix [A] we review the theory of Banach spaces of analytic functions and describe how 
it fits in our problem. In Subsection |4.4| we prove the Main Theorem and in |4.5[ we consider the 
systematic computation of higher order approximations. Besides being useful in the proof of theorems, 
the higher order approximations are the basis of efficient numerical algorithms presented in Section [6j 
In particular, in Example 6.2, we show that the method can be used to approximate slow manifolds, 
even in the presence of a singular limit. 



Some more refined analysis also leads to the study of singular limits in Section 5.1 Since the 
main tool is the implicit function theorem, we can obtain easily smooth dependence on parameters (see 



Subsection 5.2 ). 



2 General Setup 

2.1 Parameterized solutions 

We are interested in extending the theory of invariant manifolds associated to a hyperbolic fixed point 
to the more general context of difference equations. The extension of fixed point is clear. 

Definition 2.1. If 9* G M satisfies Z(8*, . . . ,9*) = 0, then we will say that 9k = 9* is a fixed point 
solution. 

The key to the generalization of the invariant manifolds from dynamical systems to difference 
equations is to observe that the invariant manifolds of a dynamical systems are just manifolds of orbits. 
This formulation makes sense in the context of difference equations. Note also that the dynamics 
restricted to the invariant manifold is semi-conjugate to the dynamics in the whole manifold (in the 
language of ergodic theory, it is a factor). Hence, we define: 

Definition 2.2 (Parameterized stable solution). Let 9* be a fixed point solution to the difference equa- 
tion ([!]). Let T> be an open disk ofM. m around the origin. We will say that a smooth function P :T> — >■ M 
is a stable parameterization of dimension m with internal dynamics h : D — > h (T>) when: 

a) p(o) = 9*. 

b) is an attracting fixed point of h. 

c) If z E V then 

Z (p (h k (zj\ , P [h k+1 {z)^ ,...,P (h k+N (z)^ = 0, (6) 

for all k > 0. 
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Remark 2.1. Notice that, in the definition above, if we let zq € T> and Ok = P(h k (zo)), then the sequence 
(9k)k>o satisfies the difference equation ([I]) and also 0^ — > 0*, as k — > oo. 

In some situations it might be useful to consider a geometric object associated to the parameteri- 



zation. As a consequence of definition 2.2, we have the following result. 

Proposition 2.3. Let 0* be a fixed point solution and P : T> C M m — > M a parameterization of 
dimension m with internal dynamics h : D — )• h(T>). If P'(0) has rank equal to the dimension ofD, 
then there exists an open disk T> Cf around the origin such that 

W = {{P(z),P (h(z)) ,...,P [h N -\z))) :zGV} (7) 

is an embedded submanifold W ^ M N of dimension m. In such case, we will say that a manifold W 
as in (m) is a local stable manifold with parameterization P. 

We note that the parameterization P and the internal dynamics h satisfy 

MP)] (z) := Z (P(z),P(h(z)), . . . , P(h N (z))) = 0, (8) 
together with the normalizations obtained by a change of origin 

P(0) = 0*, h(0)=0. (9) 

The equation ^ is the centerpiece of our analysis. Using methods of functional analysis we can 
show that, under appropriate conditions, ^ has solutions. In order to do this, we will find it convenient 
to think of Q as the equation $(P) = 0, defined on a suitable function space. 

The solutions of ^ thus produced generalize the familiar stable and strong stable manifolds but 
also include some other invariant manifolds associated to non-resonant spaces including, in some cases, 
the slow manifolds. 



2.2 An important simplification 

If L : T> — > T> is a diffeomorphism and we define P = P o L and h = L^ 1 oho L, then (P, h) solves Q 
on the domain T>. In other words, the P and h solving ^ is not unique. Nevertheless, the range of P 
and P is unique. 

One can take advantage of this lack of uniqueness of ^ to impose extra normalization on the maps 

h. Recall that, by the theorem of Sternberg [St e55l I5te58| l5Te59| . under non-resonance conditions on 

the spectrumj^] any analytic one dimensional invariant curve tangent to a contracting eigenvector has 

dynamics which is analytically conjugate to the linear one. So that, we can take h to be linear if we 

assume non- resonance or if we deal with 1-dimensional manifolds. 

lr This is automatically satisfied for the strong stable and strong unstable manifolds and includes the classical stable 
and unstable manifolds. 
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In this paper, we will not consider the resonant case. This justifies that in ([8]) we do not consider h 
as part of the unknowns, since it will be linear. In the cases considered here, it will suffice to consider 
h to be linear, we will write h(z) = Az. Hence, instead of ([8]), we will consider 

MP)} (z) := Z {P(z),P(Az), P(A N z)) = 0, (10) 

where A is a matrix that will be determined from a linearized equation near the fixed point. In fact, 
A is a matrix that depends on the roots of a polynomial. The determination of A and -P'(O) will be 
studied in Section [3| We will call $ the parameterization operator. The domain of the operator $ is a 
function space that depends on the analytic properties of Z and will be studied in detail in Section |4| 
In summary, the basic ansatz that we propose is that there are solutions to the difference equation 
that are of the form 

8 k = P (A fc z) , 

where A is determined by the linearization at the fixed point. In addition, A and P are chosen so that 
Ok — > 6* as k — > oo. The smoothness of P depends on the smoothness of Z. 

2.3 Unstable manifolds 

With the parameterization method we can also study unstable manifolds. 

Definition 2.4 (Parameterized unstable solution). Let 6* be a fixed point solution to the difference 
equation ([!]) . Let T> be an open disk of R m around the origin. We will say that a smooth function 
P : T> —> M is a unstable parameterization of dimension m with internal dynamics h : T> — >■ h (D) 
when: 

a) P(0) = 6*. 

b) is an attracting fixed point of h. 

c) If z G 2? then 

Z (P (h k + N (z)) , P (hM-^zj) , . . . , P (h k (z))) = 0, (11) 

for all k > 0. 

Remark 2.2. As in the stable case, each parameterization produces sequences that satisfy the original 
difference equation. In this case, if we let zq G T> and 6L& = P(h {zq)), then the sequence (6k)k<o 
satisfies the difference equation ([!]) and also 6k — > 0*, as k — > — oo. 

This corresponds to an extension of the original difference Z (6k, ■ ■ ■ , 6k+N) = 0, to negative values 
k < 0. Alternatively, we can write the difference equation for negative values in terms of a dual problem 

Z (&k, 6~k+i, ■ ■ ■ , 6k+N^j = Z (&k+N, 0~k+N-i, ■ ■ ■ , 6k^j = 0. 

We interpret the new variable as 6k = 6^_k- I n this way, the unstable case is reduced to the stable 
case. 
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2.4 Explicit examples 

Example 2.1. Let 77 > 0. Consider the function Z : M? — > R given by 

Z(9q, 61,62) = 6 2 + 6 ^ 2 1 • 

The resulting difference equation Z(0k,6k+i,6k+2) = [McM71] is an explicit recurrence called the 
McMillan map. It is integrable in the sense that J(x,y) = x 2 y 2 + a; 2 + y 2 — 2cosh(r])xy satisfies 
J(6k,6k+\) = J(9k+i,9k+2)- In [DRR98j . it was shown that P(z) = 2 sinh(?7)2:/(2! 2 + 1) satisfies 

Z(P(z),P(Xz),P(X 2 z)) = P(A 2 z) + P(z) - ^^jfff - 0, 

provided A = e~ ri . Therefore P is a parameterized solution with internal dynamics h(z) = Xz. Notice 
that, in the case of explicit recurrences, the parameterization can follow the manifold even if it folds 
and ceases to be a graph. 

Example 2.2. Let M = (— 00, 1). On M x M, we define the following the Lagrangian S(9o,9±) = 
\ (60 — f(6i)) 2 , where f{9) = 29/(1 — 9). The corresponding Euler-Lagrange difference equation can 
be written as Z(9k, #fc+i, ^fc+2) = where Z : M 3 — > K is the function given by 

Z(9 , 0i, 2 ) = (/(0i) - 6 ) f(9x) + (0i - f(9 2 )) . (12) 

Clearly, 9* = is a fixed point solution. 

Let P(z) = z/(l — z). Since P(z) € M, the maximal disk in which P can be defined is T> = 
(—1/2,1/2). If we let A = g; then one can check that P satisfies f(P(Xz)) = P(z), for all z G P. 
After substitution, we verify that Z(P(z), P(Xz), P(X 2 z)) = 0. Therefore, P is a stable parameterized 
solution for 0* = 0. 



However, the difference equation (12) is not a dynamical system because it can not be inverted, 
unless / is a diffeomorphism of M. In other words, we can find a parameterization solution even if the 
system is not dynamic. 

Example 2.3. Suppose that G : M d — > M d is a map, possible non-invertible, with G(0) = 0. We would 
like to parameterize the stable and unstable manifolds of 0, if they exist. For instance, we can solve the 
following pair of one-dimensional problems. 

• Stable manifold problem: if < A < 1 is an eigenvalue of DG(0), find a function P such that 
P{Xz) = G{P{z)). 

• Unstable manifold problem: if /i > 1 is an eigenvalue of DG(0), find a function P such that 
P(fiz) = G(P(z)). 

Both problems can be solved using the parameterization method for difference equations. The same 
ideas can be used for higher dimensional objects. For examples, for d = 3, two-dimensional stable and 
unstable manifolds were found in [MJL10J. As we will see, our analysis covers not only one-dimensional 
stable and unstable manifolds, but also other non-resonant manifolds. 



S 



3 Linear analysis 



Since we are interested in solutions of (10) it is natural to study the behavior of the linearization. As 
it happens in other situations, this will lead to certain choices. In subsequent sections, we will show 
that the once these choices are made (satisfying some mild conditions), then there is indeed a manifold 
which agrees with these constraints. 

Suppose that there exist differentiable P, and a matrix A solving (10). Taking derivatives of Q 
with respect to z, and evaluating at z = 0, we get that 

N 

^BiVA^O, (13) 
where Bi = diZ (9* , ... ,9*) and V = -P'(O). Our first goal is to understand conditions that allow to 



solve equation (13) which is a necessary condition for the existence of differentiable solutions for (10). 

Remark 3.1. The dimensions of the matrices are determined by the type of parametrization that we 
have. Notice that each Bi is d x d, A is m x m and V is d x m. 

Definition 3.1 (Characteristic polynomial). Let 9* a fixed point solution. If we let 

Bi = diZ(9*,...,9*), 
then the characteristic polynomial of the fixed point is defined to be: 



J"(A) := det (j^ 



If X is a root of T , then we will say that is it is an eigenvalue. The set of eigenvalues is the spectrum, 
denoted by a (9*), of the fixed point. If X G a (9*), then any vector v € C d \ {0} that satisfies 

N 

X'BiV = 

i=0 

will be called an eigenvector of X. In addition, if V is a d x m matrix with non-zero columns and 



A is a m x m matrix that satisfy the linear relation (13), then we will say that the pair (V, A) is an 
eigensolution of dimension m. 

Remark 3.2. If the difference equation is of the form ^ and is the Euler-Lagrange of a variational 
principle then the characteristic polynomial is of the form J-(X) = X Nd C(X) where C is given by 

£(A) := det j Yj xj ~ iA ij] , (14) 

y,i=o / 

and A{j = dijS (9* , . . . ,9*). Now, since Afj = Aji, we have that if A is in the spectrum, then 1/A is also 
in the spectrum. This is a well known symmetry for the spectrum of symplectic matrices. In |Ves91| 
one can find a proof that, when the Euler-Lagrange equation ^ defines a dynamical system, it becomes 
symplectic. 
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Remark 3.3. Let (V, A) be an eigensolution as above. Suppose that w is an eigenvector of A with 



eigenvalue A. Then v = Aw is an eigenvalue with eigenvalue A in the sense of definition 3.1 



Proposition 3.2. Let {y\, . . . , v m } be a linearly independent set of eigenvectors with distinct eigen- 
values Ai,...,A m . Let V be the d x m matrix V = (v\V2 ■■■ v m ) and A be the m x m matrix 



A = diag(Ai, . . . , A m ). Then (V, A) is an eigensolution and satisfies equation (13). In addition, if 



Q is an invertible matrix, V = VQ and A = Q 1 AQ then (V, A) is also an eigensolution. 

Remark 3.4. Suppose that A = /i + i^isa root of the characteristic polynomial T of a fixed point 9* . 
Let v = r + i s be an eigenvector of A. If we want to avoid the use of complex numbers, we can use A 
and v in order to find a solution with dimension m = 2. Let V be the d x 2 matrix given by V = (r s) 
and 

*-(::)■ 

It is easy to see that (V, A) is an eigensolution of dimension 2. 

Remark 3.5. We notice that J~(X) is always a polynomial of degree at most 2Nd. If the matrix Bq is 
non-singular then the degree of J-(\) is exactly Nd. If A is an eigenvalue, then there is only a finite 
number of values of n for which J-"(A n ) = 0. These considerations motivate the following. 

Definition 3.3 (Non-singularity condition). We will say that a fixed point 9* is non-singular if the 
corresponding characteristic polynomial satisfies: J-(0) ^ 0. 

Definition 3.4 (Non-resonance condition). We will say that an eigenvalue A G a (9*) is non-resonant 
if F{\ n ) ^ 0, for alln>2. 

More generally, we will consider non-resonant sets of eigenvalues. In what follows, we will be using 
the multi-index notation where, as usual, if z = {z\, . . . , z m ) £ C m and a = (at, . . . , a m ) G is a 
multi-index then z a = z" 1 z^ 2 ■ ■ ■ z^ 1 . 

Definition 3.5. We will say that X = (Ai, . . . , A m ) is a non-resonant vector of eigenvalues if, for all 
multi-indices a G Z™, 

a) T(X a ) = if\a\ = I, 

b) T(X a ) ^ if\a\ > 1. 

If X = (Ai, . . . , A m ) is a non-resonant vector of eigenvalues then we will also say that the set {Ai, . . . , A m } 
is non-resonant. 

Definition 3.6 (Hyperbolicity). Let 9* be a non- singular fixed point solution of an analytic difference 
equation Z and suppose that T is its characteristic polynomial. We will say that 9* is hyperbolic if 
none of the eigenvalues of J- are on the unit circle. Similarly, we will say that a vector of eigenvalues 
X = (Ai, . . . , A m ) is stable if |Aj| < 1, for all i = 1, . . . ,rn. 
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Remark 3.6. Let 9* be a non-singular fixed point solution and cr(6*) its spectrum. In other words, all 
elements of <j{9*) are non-zero. Notice that, even if the condition of non-resonance involves infinitely 
many conditions, for stable sets all except a finite number of them are automatic. Suppose that 
A = (Ai, . . . , A m ) is a stable vector of eigenvalues. If n £ N is such that 

( Y 

( max I AJ < min |A|, 
\\e\ J Aea(0*) 

then we have that X a cannot be in the spectrum if |a| > n. So that there are only a finite number of 
conditions to check and the set of non-resonant A is an open-dense, full measure set among the stable 
ones. 

All this analysis shows that there are obstructions to the computation of invariant manifolds that 
appear using just the linear approximation. These obstructions are a generalization of the observation, 
in the dynamical systems case, that the only invariant manifolds have to have tangent spaces that are 
invariant under the linearization. 

The goal of the rest of the paper is to show that if we choose a subset of the spectrum and an 
invariant subspace, which also satisfies the non-resonance conditions, then there is a solution to the 
parameterization problem. The analysis will also show that the non-resonance conditions are needed 
to obtain a general result. 



4 Existence and analyticity of solutions 

As indicated before, we will first obtain a formal approximation and then use an implicit function 



theorem. The proof of Theorem 4.1 is an application of the implicit function theorem in a Banach space 



of analytic functions (see Subsection 4.4). This will lead rather quickly to the analytic dependence on 



parameters, and the possibility of getting approximate solutions (see Subsection 4.5). 



4.1 Statement of the Main Theorem 

Theorem 4.1. Let Z be an analytic difference equation function with a fixed point at and such 
that Bq = doZ(0) is a non-singular matrix. Let A = (Ai, . . . , A m ) be a stable non-resonant vector of 
eigenvalues, with corresponding eigenvectors v\, . . . ,v m . Let (V, A) be the corresponding eigensolution 



of the linearized problem (13) with A = diag(Ai, • • • , A m ) and V = (v\ ■ ■ ■ v m ). 

Then, there exist an analytic function P such that its derivative at z = is -P'(O) = V and satisfies 
Q and (10). The solution is unique among the solutions of the equation. 



Remark 4.1. The radius of convergence of the resulting analytic function is not specified. In the proof 



of 4.1, we will consider an equivalent result, in which the radius of convergence is 1, but the lengths of 
the vectors v\, . . . , v m are modified. 

Remark 4.2. If Z is the Euler-Lagrange equation ^ that corresponds to a generating function S 1 , then 
-Bo = do,NS(fy an d hence the non-singularity condition is det(9o,jv£(0)) ^ 0. 
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In Section 5.1 we will weaken the assumption det(-Bo) 7^ 0, which is tantamount to a local version 
of the Legendre condition. This includes, in particular, the extended Frenkel Kontorova models with 
singularities. 



4.2 Analyticity of the parameterization operator 

In appendix |Aj we give a general definition of analyticity and introduce spaces of analytic functions. 
In particular, given p > we define 

( 



A p (C e 



/ : C £ (p) -> C c 



/(*) 



EE 

k=0 \a\=k 



< oo 



where 



E m« 

k=0 \\a\=k 



P 



|| • Hoc is the uniform norm (45), and C (p) = {z £ C : ||^||oo 5; p}- 

Let d, N £ N and £ = d(N + 1). We will consider C e ~ (C d ) JV+1 
the function that defines the difference equation satisfies Z G A p (C 
analytic near the origin. 



As an standing assumption, 
d ) for some p > 0; this is, Z is 



In this section we will show that the operator <I>, defined in (10), is an analytic operator defined on 
spaces of analytic functions in C m . Let A be an m x m matrix such that, if ||z||oo < 1 then HA^loo < 1. 



We define the linear function Q : A\ 



•"m /nd 



Q{P) 



P,PoA, 



C e ) by 
,PoA N ) 



Besides being linear, the functional Q has the property that Q [A{ (C m ,C d )) C A{ (C m ,C^) , for all 
p > 0. Let 

X = A 1 (C m ,C d ) (15) 

and define the open set 

U = {feX:\\f\\ 1 <p}. (16) 

For a correct application of the implicit function theorem, it suffices $ to be C . However, we next 
show if is analytic. This will give analytic dependence on parameters of the invariant manifolds. 

Notice that Q (U) C A^ (<C m , C^) . Since Z £ A p (<C^, C d ) , we can write the parameterization operator 
as a composition of analytic functions: $ = Cz°Q, where Cz is the composition operator Cz{g) = Zog. 
Consequently, using Lemma A. 4 we get the following: 



Proposition 4.2. Let U C X as in (15) and (16). Then the operator : U — > X is analytic. 
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4.3 Frechet derivative of the parameterization operator 



Consider the parameterization operator $ defined in (10), in which A is a diagonal matrix. As we have 



seen, this operator is a function : IA — > X, where X and IA are defined in (15) and (16) respectively. 
We have the following. 

Lemma 4.3. Let A = (Ai, . . . , A m ) be a vector such that ||A||oo < 1- Let <J> :IA — >■ X be the parameteri- 
zation operator defined in (10) where A = diag(Ai, . . . , A m ) and define 

N 



T(A) :=^A^, 



(17) 



i=0 



where B, L = diZ(0, . . . , 0). Let p G X be of the form <p(z) = ^ V] z a ip a . Then D<fr(0)ip G X and 

k=0 \a\=k 



k=0 \a\=k 



(18) 



Proof. We notice the conditions on A imply that ip o A* el for alii = 0, . . . , N . We also have that the 
Frechet derivative of satisfies 



N 



(19) 



8=0 



for all Hzlloo < 1. This implies that D$(0)p G X. Clearly, (A*z) a = z a (A a )' and therefore 

oo 

<p(a*z) = E z " ( Aa )^- 

fc=0 |a|=fc 



Combining the last equation with (17) and (19), we get (18). 



□ 



Let H be the Banach subspace of analytic functions in the unit disk, that vanish at the origin along 
with their first derivatives. 



H = < P> el 



^)=EE^ 11^111 = EE ii^ 



< oo 



(20) 



k=2 \a\=k 



k=2 \ a \=k 



Remark 4.3. In the notation of the appendix |a| this subspace is just H = {P G Ai(C m ,C d ) : P(0) 
0,P'(0) = 0}. 



Clearly, (18) implies that H is invariant under _D<J?(0). In addition, we get the following result. 

Lemma 4.4. Let A = diag(Ai, . . . , \ m ), where A = (Ai,...,A m ) is a non-resonant stable vector of 
eigenvalues. If Bq is a non-singular matrix, then 

a) D<&(0) is invertible in H , with bounded inverse. 

b) If cp E H is such that p(z) = ^2^L 2 Yl\a\=k za ¥a and n = DQ(0)~ 1 p> can be written as r/(z) = 
Y1T=2 Yl\a\=k za Va, then rj a = T(\ a )~ ip a , for every multi-index a such that \a\ > 2. 
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Proof. Clearly, as k — > oo we get that T(X k ) — > Bq, a matrix that is invertible. This implies that there 
exists a constant Co > and a radius < 5 < 1 such that, if |A| < S then ^(A)^ 1 !! < Co- In particular, 
if |A a | < 5 then 1 1 TT-X" ) 1 1 1 < C . 

We know that, since A is non-resonant, the matrices T(A Q ) are invertible, whenever \a\ > 2. Since 
A is stable there is only a finite number of elements in the set 

{X a G C : |a| > 2,|A a | > 5}. 

Let C > be a constant such that C > Co and 

C > max{||T(A Q )~ 1 || : \a\ > 2, |A a | > 6}. 

This implies that ||r(A Q ) _1 || < C, for all \a\ > 2. 

From ( [l8| ) we get that L><3?(0) is injective. Let <~p G iJ with </j(z) = £|a|=fc zC Va- Define 



r /(- z ) = Sfc^=2 S|a|=fc 'Ja 2 '*! with 7/qj = T(\ a ) (fa, for every multi-index |a| > 2. From (18), we have 
that D<fr(0)rj = <p and therefore Z?$(0) is invertible in H. In addition, 

00 00 

||i>*(o)-V||i = Wvh = EE l|r(A Q )-Valloo < ii^iu = cii^Hi. 

k=2 \a\=k k=2\a\=k 

We conclude that D®(0) is invertible in H and j | ^><E> (0) 1 1 1 <C. □ 
4.4 Proof of Theorem 14. 1L 



Let H as in (20). When we consider the linear part V of the parameterization P, it is convenient 
to choose the scale sufficiently small. Choosing this scale is tantamount to choosing the radius of 
convergence of the solution P. Let r = (t%, . . . , r m ), where T\, . . . , r m represent scaling factors for the 
columns of the linear part. We will use the notation r • V = Vdiag(ri, . . . , r, 



Since any possible solution of ( 10 ) has to match the lower order terms found, it is natural to consider 
a new decomposition P = r • V + P > where P > is an analytic function which vanishes to order 2 and 
depends on the size of the scale r. Because of the change of variables, we can seek for P > among 
analytic functions of radius 1 that vanish to first order, i.e. P > G H. 



Hence we write the equation ( 10 ) as 

*(r,P > ) = $(r- V + P>) = 0. (21) 

It is important to note, since V is known, that we only need to find the appropriate scale r, and the 
function P > . Furthermore, since P > vanishes to order 1, \P(r, P > ) also vanishes to order 1. In other 
words, we can choose H to be the codomain of \&. In addition, the coefficients of P > are small if r is 
small. 

We notice that there exists an open subset V C M m x H defined by the property that if (r, P > ) G V 
then t • V + P > G U, the domain of In this way, we can restrict the domain of $ and consider 
: V — > H. It is clear that this is a neighborhood of the origin (0, 0) in which the operator \I' is defined 
and is analytic. 
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Taking the derivative with respect to the second variable, it is clear that .02^(0,0) = D$>(0) and 
therefore, by Lemma 4.4 we have that the operator 



D 2 *(0,0) : H -»• H 

is invertible with bounded inverse. The implicit function theorem in Banach spaces [HS91J implies that 
there exists 5 > and a function r i— > P^~ such that if \t\ < 8, then ^(r, Pf*) = and the solution is 
unique if we require ||-P > ||i < 5. 



Fix a solution of (21) of the form (t, P>) such that r has positive entries and ^(r, P^) = 0. Then 
P(z) = t -V z + P^ (z) is a solution of the parameterization problem. By construction, P is an analytic 
function with radius of convergence 1. If we want to have V as the linear part of the solution, we modify 
P in the following way. Let P = P o diag(ri, . . . , r m ) _1 or 

P{z) = Vz + P>(diag(n, . . • ,r m )^ 1 z). 



Using a linear change of variables as in Section 2.2 we conclude that P is also a solution of the 



parameterization problem. However, the radius of convergence is not longer 1 but depends on the choice 
of t. If we let r = minjri, . . . ,r m }, then r > 0, and ||z||oo < r implies that || diag(ri, . . . ,T m )~ 1 z\\ oc < 1. 
We conclude that P E A r {£, m ,C d ) and has radius of convergence equal to r. □ 

4.5 Formal approximations to higher order 

Once -P'(O) = V and A are chosen, the solution P of the parameterization problem can be approximated 
with the first terms of the power series. Due to analyticity, we can write the solution P G Ap(C m ,C d ) 



of the parameterization problem as a sum of homogeneous polynomials like in ( 50 ) : 



l=\ \a\=l 

where P a G C d , and has P radius of convergence p. 

For each multi- index a £ , we will use the notation [•]„ for the coefficient vector of the term that 
corresponds to z a . Clearly, if / is an analytic function at the origin, then this coefficient can be written 
as [f] a = d a f(0)/al. In the case of P above, we get that [P] a = Pa- 

For each n G N, let P- n be the polynomial 

n 

1=1 \a\=l 

It is clear that \\P — P- n ||p — > as n — > oo. Notice also that, for any \a\ < n, 

[$(P<-)] Q = [$(P)] a = 0, (22) 

We will describe how to construct the polynomials P- n recursively, provided that the eigenvalues are 
non-resonant. 
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4.3 



implies that 



A simple computation shows that [<J> (P- 1 )] a = for \a\ = 1. Let Af(P) = P$(0)P - <&(P) 
Then, it turns out that [M {P- n )] a = [N { p - n+1 )) a > for a11 \a\ = n + 1. Lemma 
[D$(0)P] a = T(\ a )P a , for all \a\ > 1. Therefore, we conclude that 

0= [$ (P^" +1 )] a = [P$(0)P^ n+1 +A/'(P^ n+1 )] Q 
= [ J D$(o)P^™+ 1 ] Q +[AA(P^)] a 
= T(X a )P a +[M(P^ n )] a , 

From this we find the expression 

p a = T(\ a )~ 1 [at (p* n )] a = r(A a )- 1 [$ (P* ")] a , 

for all \a\ = n + 1. The polynomial _P- n+1 can be found from P- n and recursion (23) 



(23) 



It is important to note that the choice of P'(0) determines the tangent space to the manifold. Hence 
V is determined once we choose this space. On the other hand, P'(0) is determined only up to the size 
of its columns. These multiples will not be too crucial for the mathematical analysis, but it will be 
important in the numerical calculations in Section [6j 

Lemma 4.5. With the notations above, assume that A = (Ai, . . . ,X m ) is a stable non-resonant vector 
of eigenvalues, and that (13) is satisfied with A = diag(Ai, A m ) and for some V = P 1 (0) of maximal 



rank. Then, for every a > 2, we can find a unique P a such that (22) holds. Furthermore, we can 



make all P a arbitrarily small by making the columns o/P'(0) sufficiently small. 

Remark 4.4. The assumption that the matrix A is diagonal can be eliminated. Following the discussion 



in Section 2.2 it suffices that A is diagonalizable. 



Remark 4.5. As we will see in Section [6] the proof in this section can be turned into an efficient algorithm 
using the methods of "automatic differentiation" [JZ05J INeilOl |BCH + Q6 which allow a fast evaluation 
of the coefficients P a , specially in the case that the manifolds are 1-dimensional. More details, including 
an implementation in examples, are given int Section [6} 

Remark 4.6. Notice that the main theorem is proved by a contraction mapping theorem. The formal 
solutions are indeed an approximate solution. Indeed, in practical problems - see Section [6] - it is 
possible to produce solutions that have an error comparable to round off. These bounds can be proved 
rigorously using interval arithmetic. 

Given some bounds on the contraction properties of the operator, one concludes bounds on the 
distance between the approximate solution and the true solution. Hence, the proof presented here gives 
a strategy to lead to computer assisted proofs. 



5 Singular case and dependence on parameters 
5.1 Singular case 

In this section, we show how the results can be extended to the case in which Bq is singular. The key 



will be to generalize Lemma 4.4 This can be done by estimating the singularity of the matrix T(A) 
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that was defined in (17). 

Example 5.1. Consider the Lagrangian function S : M 2 x M 2 — > K given by 

S(0a,9i) = 

Let Z be the difference equation that arises from the Euler-Lagrange equation ([3]). The point 8* = (0, 0) 




gives a fixed point solution. As before, we define T(A) as in equation (17). Using definition ( |14[ ) in 
remark 3.2 it is possible to verify that the characteristic polynomial is of the form J-(X) = A(— 2A + 
1)(A — 2). Therefore, 8* is singular. We notice that the degree is strictly less than the maximum Nd = 4, 
J"(0) = and A divides J"(A). 

In order to deal with the singular case, we have the following result, that will lead to a generalization 
of Lemma 14.41 



Lemma 5.1. Let 8* be a singular fixed point solution with characteristic polynomial T . Let e be the 
greatest integer e £ Z + such that A 6 divides J~(X), i.e., the polynomial J- is of the form J-(X) = X e g(X), 



where g is a polynomial such that g(0) ^ 0. Let T be as in (17) and X = (Ai,...,A m ) be a stable 
non-resonant vector of eigenvalues none of which is zero. Then there exists a constant C > such that 



II^a")- 1 )! < c{\ a y 

for all multi-indices \a\ > 2. 

Proof. The inverse T(A) _1 is a rational function of the form 



(24) 



T(A)- 1 



1 



:Q(X), 



X e g(X) 

where Q is a polynomial matrix and g(0) 7^ 0. Then A e T(A) _1 is also a rational function, but the limit 
lim^o A e T(A) _1 exists. 

As in the proof of Lemma |4.4[ we can argue that since A is non-resonant and stable, the matrices 
T(X a ) are invertible and (A a ) e T(A a ) -1 are uniformly bounded for all multi-indices \a\ > 2. Therefore, 
there exists a constant C > such that ||(A a ) e T(A Q )~ 1 || < C. □ 

Lemma 5.1 tells us that the derivative D$(0) _1 is an operator, but it might not be well defined for 
all the elements of H. We will introduce a new Banach space. For each /j, > 0, let D(fi) = C m (fj,) = 
{z e C m : \\z 
define 



< fi} be the complex disk around the origin of complex radius 11. Using D(fi), we 



H(jm) :-- 



P> : D(fi) -»■ C a 



k=2\a\=k k=2 \a\=k 



Remark 5.1. Notice that, in the notation of the appendix, 

H(ji) = {P£ ^(C m , C d ) : P(0) = 0, P'(0) = 0}. 
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If fii < fi2 then D(m) C D{^%) and H(fi2) can be regarded as a subspace of H(jx\) through the 
standard inclusion H(fj,2) > H{l^i) given by 

/!%)■ 

In particular, we have that H = H(l) and if < \i < 1 then we have the inclusion H ^ H(fi). In that 
case, we define an operator A M : H(fi) — >• H by 

A ^M(^) = </>(/^)- 

Clearly, is always bounded. Functions in H(fi) gain analyticity through A^, and the inverse A^ 1 is 



also bounded but destroys analyticity. Following the proof of Lemma 4.4 as a corollary of Lemma 5.1 
we have the following result. 

Corollary 5.2. Let 9* be a fixed point solution and A = (Ai, . . . , A m ) be a stable non-resonant vector 
of eigenvalues. Let e be the greatest integer e € such that A 6 divides J-(X) . 

Suppose that < [i < min{|Ai| e , . . . , |A m | e }. Then D^O)^ 1 is a bounded linear operator D<I>(0) _1 : 
H — > H(ijl). In addition, the composition A^D^O) -1 : H — )• H is bounded and 

llA^tO)- 1 !! < C, (25) 



where C is any constant that satisfies (24) in Lemma 5.1 



5.2 Dependence on parameters 

Suppose that the difference equation depends smoothly on q parameters that, for simplicity belong to 
an open set E C around the origin. A special difficulty arises when there is a parameter in which 
the equation becomes singular. We would like to regularize the singular limit. 

Example 5.2. Consider the Lagrangian S defined on M 3 and given by 

S(0 O , 01, 02) = | (02 - O ) 2 + #0) 2 + W(0 O ), 

where e is a small parameter. From this, we get the Euler-Lagrange equation expressed as (|3|. If 
W'(0) = 0, then the point 0* = (0,0) is a fixed point solution. The characteristic polynomial for this 
point is 

J-(A) = eA 4 + A 3 - (2s + 2 + W"(0)) A 2 + A + e. 



According to Definition 3.3 the fixed point 0* = (0, 0) is non-singular if J~(0) ^ and this happens if 
and only if e ^ 0. 

As before, let X = Ai(C m ,C d ). As an assumption, suppose that there is an open set U C X such 
that the equation : £ x U — > C d is analytically defined. Suppose that 0* = G U is a fixed point 
solution for all parameters. On £ x U, we define the nonlinearity at 0* as 

M(e,P) = <Z>(e,P)-D 2 <Z>(e,9*)P, 

for all (s, P) € £ x U. 
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Remark 5.2. In general, for each value of e G £, the corresponding characteristic polynomial J- £ (\) is 
of the form ^(A) = X e ^g £ (X), where e(e) is an integer that depends on e and g £ (0) 7^ 0. The function 
e(e) does not need to be continuous and this constitutes a potential difficulty. However, the theorem 
below only requires that one can find a constant fi for which the matrices /i' Q 'T(A(0) a ) _1 are uniformly 
bounded and the nonlinearity N vanishes at higher order. 

The solution of the stable manifold problem is a local issue, so we can consider IA and £ as small 
neighborhoods that not necessarily cover the largest possible domain for <£. Now we can prove a more 
general theorem that includes parameters. 

Theorem 5.3. Let <£, 8* , N ', £ andlA as above. Suppose that there exists analytic functions Aj : £ — > C, 
Vi : £ — > C d , for i = 1, . . . ,m and constants C, /x > such that the following conditions are satisfied, 
for all e G £ . 

a) Each Aj(e) is a non-resonant eigenvalue with eigenvector Vi{e) . 

b) A(e) = (Ai(e), . . . , A m (e)) is a stable non-resonant vector of eigenvalues. 

c) ||T(A(0) a ) _1 || < C/i~' a ', for all multi-indices such that \a\ > 2. 

d) There exists an open neighborhood of the origin Uq C Li such that the operator lZ(e, P) = 
M (e, A^ 1 P) can be defined as a function 1Z : £ x Uq — > X . 

For each e G £, let V(e) = (vi(e) ■ ■ ■ v m (e)) and A(e) = diag(Ai(e), • • • ,A m (e)). Then, there exist a 
function P e {z), analytic in z and e such that P £ (0) = 0, its derivative at z = is d z P' £ (Q) = V(e) and 

Z £ (P £ (z),P £ (A(e)z), . . . , P £ (A(e) N z)) = 0, 

for all z in a neighborhood of the origin. The solution is unique among the solutions of the equation. 

Proof. For simplicity, let the fixed point solution be 0* = 0. We will use the notation <3? e = <l>(e, •). 

Let H as in ( 20 ) . The main step of the proof is to show that there exists 5 > and an analytic 
function (r, e) h-> Q> e defined for all HrHoo + ||e||oo < <5 such that 

<5>(e,r-V(e) + Q> £ ) =0, 

with Q> £ G H and the solution is unique provided ||Q^ e ||i < 6. 

Now, since G Uq, there exists an open set V C W 71 x £ x H that contains (0,0,0) such that if 
(r,e,Q > ) G V then r • (fiV(e)) + Q y G Uq. From equation ([25]), we get that D$(0) _1 is a bounded 
linear operator L><I>(0) _1 : H — > H(/i) and A M D<I>o(0) -1 : H — >• H is uniformly bounded. 

Let TZ(e, P) = N (e, A^ 1 P) . Notice that, if (r, e,Q > ) G V then TZ(e, r • (jJtV(e)) + Q > ) G H. Let 
^ : V — > H be the operator defined by 

*(t, e, Q>) = J D$o(0)A- 1 Q > + 7J(e, r • (pV{e)) + Q>). 
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We notice that * (0,0,0) = and D 3 *(0, 0, 0) = D$ (0)A~ 1 , that by construction is invertible with 
bounded inverse. Using the implicit function theorem of Banach spaces, we can find 5 > and a 
function (t, e) \-t Q> e £ H defined for HtHoo + ||e||oo < S such that \&(t, e, Q^ e ) = and the solution is 
unique if ||Q^ e || < 5. For each such (r, e), define P> e = A~ Q> e . This implies that 

D<f> £ (0)P> e + 7V(e, r • V(e) + P> £ ) = 0. 

Since D$> £ (0)V(s) = 0, we have that P T:£ = r • V(e) + P> £ is a solution and satisfies 

<D e (P T , e ) = <fr e (r • F(e) + P> ) = 0. 



The proof is finished as in the proof of Theorem 4.1 with a change of variables. Suppose that P^ e is a 
solution such that the vector r = (7*1, ... , r m ) satisfies || 7~ I J 00 ^ t ciiicl its entries are positive. Let 

P £ = V(e) + P> e o diag(n, . . . , g- 1 . 

Then P e (0) = 0, 5 2 P^(0) = U(e), and has a radius of convergence r, where r is given by 

r = min{n, . . . ,r m }. 

□ 



6 Examples of numerical algorithms 

In this section we describe efficient numerical methods to compute the parameterizations P described 
in the previous sections. We will present algorithms for 

A. Standard map model with several harmonics. 

B. Frenkel-Kontorova models with long range interactions^] 

C. Heisenberg XY models. 

D. Invariant manifolds in Froeschle maps 

For simplicity, the invariant manifolds are 1— dimensional. System A is a twist map, B contains a 
singular limit, C does not define a map. In D is a 4— dimensional and we study both strong stable and 
slow invariant manifolds. 

6.1 Standard map model with K— harmonics 

Let Ci, . . . , Ck be given numbers. Consider the Lagrangian S : M 2 — > M. given by S(9o,6\) = ^(#0 — 
#i) 2 + W(9 ), with 

K 
j=Q 3 

2 We will have the C code available in the web. 
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The the corresponding Euler-Lagrange difference equation can be written as Z(6q,9i,02) = 0, where 
Z : M 3 -> R is the function given by Z(9 , 1 ,9 2 ) = 9 2 - 20 1 + 9 - W'{9i). 

Many authors have treated the original (K = 1) standard map, also known as Chirikov model. The 
bi-harmonic model (K = 2) was studied by [BM94 } lLC06| . The parameterization problem to solve is 



A' 



P(X 2 z) - 2P(Xz) + P(z) - ^ Cj sm(jP(Xz)) = 0, 



(26) 



j=0 



where A solves P(A) = A 2 + ( — 2 + Ylf=o 3@j ) A + 1 = 0, and | A| < 1. Sometimes, it is more convenient 



to write ( 26 ) as 



A 



P(Xz) - 2P(z) + P(A _1 «) - c i sin(jP(z)) = 0. 

3=0 



(27) 



We solve (27) by equating coefficients of like terms. The orders and 1 are special. Clearly, we can 



take Pq = 0. This corresponds to choosing the fixed point 9* = 0. Equating terms of first order in (27), 
we obtain: 

A 



A + A- 1 -2 + ^jC j ] P x 

j=0 



0. (28) 
Since we choose A so that the term in parenthesis vanishes, we obtain that Pi is arbitrary. Once we 



have chosen A so that it solves the quadratic equation, any Pi will lead to a solution of (28). 

Any choice of Pi is equivalent from the mathematical point of view as they correspond to the choice 
of scale of the parameterization. However, from the numerical point of view it is convenient to choose 
Pi in such a way that the subsequent coefficients have comparable sizes so that the round of error is 
minimized. In practice, to find a good choice of Pi we perform a trial run of low order which gives an 
idea of the exponential growth or (decay) of the coefficients P n and then fix Pi so that the coefficients 
P n neither grow nor decay too much. 

The equation for A is quadratic. The product of its roots is 1 so, when 



A 
3=1 



> 2, 



we get two roots Ai and A2 of the characteristic polynomial T such that |Ai| < 1, | A2 1 > 1. We choose 
the stable eigenvalue Ai. 



Since X n + X~ n - 2 + Y^=i 3 Cj + (A n is not a root of P), we get 




A 



^Cj sin (jP<("-D 



(29) 



Note that the right hand side can be evaluated if we know P-( n_1 ) and hence we can recursively 
compute P n . In each step, the coefficients can be found using algorithms explained below which will be 
also used in other sections. 
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6.2 Efficient evaluation of trigonometric functions 

Given a series, P(z) = Y^^=o PnZ n , we often want to compute the power series expansions of sin(P(z)) 
and cos(P(z)). The following algorithm is taken from |Knu97j . Denote S(z) = sin(P(z)) and C(z) = 
cos(P(z)). Then, we have 

S'(z) = C{z)P\z), C'(z) = -S(z)P'(z). (30) 

Suppose that we can write these functions as S(z) = Yl^=o $nZ n and C(z) = Yl^=o C n z n ■ For each 
n € N, we will denote S^ n (z) = Y2=o S kZ k , C^ n {z) = Yl=o C kZ k and P^ n (z) = Yl=o p kZ k - Also, 
[•] n will represent the coefficient of order n of an analytic function. 



Equating terms of order n in (30), we obtain 



(n + l)5 n+ i = Cn -^ + 1 ) P i+ 1 ' ( 31 ) 

3=0 
n 

(n + l)C n+ i =-J2 s n-jU + l)Pj+i- 

j=0 



The recursion (31) allows to compute the pair of coefficients S n+ i, C n +\ provided that we know the 
coefficients So,...,S n and Co, . . . , C n . We note that, obviously So, Co are straightforward to compute. 
From this, we also make the obvious observation that 

S n+ i = CoPn+i + ^ [cos (P^ n ) Q- n ] n+1 , (32) 

C n+1 = -S P n+1 - [sin Q^] n+ i > 

where Q^ n {z) = Y2=o( k + ^)Pk+iz k - In particular, if P = 0, then we conclude that 

[S] n+1 = P n+1 +[sm(P^)] n+1 , 
[C] n+1 =[cos(P^)] n+1 . 

This recursion allows us to get the expansion to order n of sin(P(z)) and cos(P(z)), given the 
expansion of P to order n. Furthermore, we observe that if we change P n -the coefficient of order n of 
P— this only affects the coefficients of sin(P(z)) and cos(P(z)) of order n or higher. 

The practical arrangement of the calculation of the coefficients in the standard map with K harmon- 
ics is to keep different polynomials Sf n and Cf n that correspond to the series expansions of sin(^P) 
and cos(£P) up to order n. If the polynomial P is computed to order n — 1 and the Sf n and Cf n 



corresponding to P- ( n 1 ) are computed up to order n, we can compute the coefficient P n using (29) 



Then, we can compute the corresponding Sf ^ n+1 ^ and Cf ( - n+1 ^ up to order n + 1 using (31 ). We note 
that similar algorithms can be deduced for e p ( z \ log P(z), P{z)^ or indeed the composition of P with 
any function that solves a simple differential equation. 
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6.3 Frenkel-Kontorova model with extended interactions 
6.3.1 Set up 

Consider the Frenkel-Kontorova model with long range interactions. A particle interacts not only 
with its nearest neighbors, but with other neighbors that are far away. Let N > 2. We consider the 
Lagrangian function S : l^ 1 — > M given by 



1 N 

S(6 , ...,e N ) = -Y,7L(0L- ) 2 + W(9 ). 



L=l 

The corresponding Euler-Lagrange equilibrium equations are expressions of 2N + 1 variables, that 
in this case have the form: 

TV 



(o k+ L -20 k + 0k- L ) - w\e k ) = o. 



L=l 



These equations represent a difference equation of order 2N . 

Suppose that ^'(0) = 0. Then the system has a fixed point solution at and the corresponding 
characteristic polynomial is of the form ^(A) = \ N C{\), where 



TV 



£(A) = 7L (A L - 2 + \~ L ) - W"(0). 



L=l 

In addition, the one-dimensional parameterization equations of the point can be written as 

N 

Y 1L {P(\ N+L z) - 2 P(X N z) + P(X N - L z)) - W'(P(X N z)) = 0, 

L=l 

where A is a non-resonant stable root of the characteristic function C. 

Remark 6.1. We can simplify the characteristic polynomial above. Notice that, if we let oj = (A+A _1 )/2, 
then 

^ = 7l(w), 

where 72 is the L— th Tchebychev polynomial. Let r(u>) be the polynomial of degree iV given by 

N 1 
r( W ) = ^ 7i (T L M-l)--r(0). (33) 

L=l 

Then, characteristic polynomial -^(A), can be written as J"(A) = 2\ N r ((A + A _1 )/2) . In addition, 
has no zeroes on the unit circle if and only if r(u) has no roots on the segment [—1, 1] C C. For 
each root u of r, we get a pair of eigenvalues. If uj is real and \uj\ > 1, then these eigenvalues are a pair 
of real numbers 

X s ' u = oj± V^ 2 - 1 

that satisfy < |A S | < 1 < |A U |. 
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6.3.2 Singular limit and slow manifolds 

In many situations the long-range interactions of the particles in the model are small. We could ask 
the question of what happens in the limit. It turns out that the system becomes singular and the 
usual dynamical systems approach fails to be useful. However, certain stable manifolds persist, as in 
Theorem |5.3| We illustrate this difficulty with an example. 



Example 6.1. Consider a Frenkel-Kontorova equation with 71 = 1 and 72 = e. In this case, the 

r(u) = e(2oo 2 - 2) + 00 - /3, 



auxiliary polynomial in ( 33 ) is 



where /3 = 1 + \W"(G). 

Solving for 00(e) , we get that 



1 ± ^/T+8eJP + 2ej ' 

If e — > then we have a singular limit. We notice that, as e — > 0, the two roots of the polynomial 
have two different limits oo + (e) — > (3 and co~ (e) — > 00. In terms of the stability of the fixed point, the 
limit oo~ (e) — > 00 corresponds to a pair of eigenvalues A s , A" that are very hyperbolic in the sense that 
A S A" = 1 and X s -> and A" -> 00 as e -> 0. 

Fortunately, the other pair of eigenvalues can be continued through the singularity e = 0. This family 
is smooth and will be denoted by X s '(e), X u (e). They satisfy < A s (e) < 1 < X u (e), X s (e)X u (e) = 1 and 
A s (0) + A u (0) = 2/3. 

Remark 6.2. In general, if we let f3 = 1 + W" (0)/ (2ji), then there is a family of roots oj of r(oj) such 
that co — > (3 as (72, . . . , 7jv) - > 0. It follows that, if \j3\ > 1 and the coefficients 72, • • • ,7tv are small 
enough, then the fixed point 6* = is hyperbolic. This occurs, for instance, when is a minimum of 
the potential W, 71 > 0, and the long range interactions are weak. 

We are interested in the persistence of slow manifolds in the Frenkel-Kontorova model with long- 
range interactions. We can consider that the interactions are small. Suppose that we have N long-range 
interactions represented by small coefficients 72(e), ... , 7jv(e) that depend analytically on the parameter 
e. Assume that 72(0) = • • • = 7at(0) = 0, j' N (0) 7^ and, without loss of generality, that 71(e) = 1. 

For each e, the characteristic polynomial T e of the fixed point 9* = is of degree at most Nd. From 
the implicit function theorem, we can argue that there exists a number eo^O and a smooth function 
00(e) such that, if |e| < £0 then r e (uj(e)) = 0, where 00 (0) = {3 and 

N 1 

r B (w) = Y,lL{e){T L {u) - 1) - -W"(0). 

L=l 

All the other roots of r £ (uo) diverge as e — > 00. For each non-singular root oo(e) of r £ , we get the 
following stable eigenvalue 

X s (e) = 00(e) - ^oo(e) 2 - 1. 
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The family A s (e) can be continued through the singularity e = and, for all |e| < £q, it corresponds to 
a slow manifold, i.e. the invariant manifold with the largest stable eigenvalue. In fact, X s (e) is analytic 
near e = 0. 

For each |e| < £q, there exists a non-negative integer e(e) such that -F e (A) = X e( - £ ^g £ (X). By con- 



struction, e(e) = if e and e(e) = N — 2 if e = 0. Using the notation of Theorem 5.3, we have that 
the nonlinearity of the parameterization operator at 9* is precisely 

[Af(e, P)](z) = W"(0)P (\ s (e) N z) - W (P(\ s (£) N z)) . 

Let \x be a number such that 

A^O)* < ii < A^O)^ 2 . 

Restricting eq further, it is possible to assume that fi also satisfies X s (e) N < ji < X s (e) N ~ 2 , for all 



M < £o- Using Lemma 5.1, we conclude that there exists a constant C such that the matrix norm 
satisfies 

||A s (0)( iV - 2 ) fe r(A s (0) fc )- 1 || < C, 
for all k > 2. Furthermore, we have that 

\\n k T(\ a (0))- 1 \\ < Cn k X s (0)-^ N -^ k < C, 

for all k > 2. Now, the condition on the nonlinearity is that TZ (e, P) = M (e, A" 1 / 5 ) € X, for all P G X 
in a neighborhood of the origin. However, the nonlinearity satisfies 

[71 (s, P)} (z) = [AT (e, A^P)] (z) = W"(0)P X s (sf z) - W (P(^X s (s) N z)) . 

Since |^~ 1 A s (e) iV | < 1, we conclude that there exists an open neighborhood IA$ of in which the 
nonlinearity can be defined as an operator. From these considerations we conclude that the Theorem 



5.3 applies. Therefore, there exists a family of analytic solutions P £ of the parameterization problem. 



This is illustrated in the numerical example that follows. 
6.3.3 Some numerics 

If we consider the K— harmonic potential W(9) = —6 Ylj=i ~f" cos (j^)i then the equilibrium equations 
are 

N K 

lUh+L ~ 29 k + 9 k _ L ) + 6^Cj sm(j9 k ) = 0. (34) 

L=l j=l 

The parameterization equations of a 1— dimensional stable manifold can be written as 

N K 

[*(P)] (z) = £ lL {P(X L z) - 2P(z) + P{X- L z)) + Sj2 C i sm(jP(z)) = 0, (35) 

L=l j=l 

where A is a non-resonant stable eigenvalue. By symmetry, if P parameterizes a stable manifold it also 
parameterizes an unstable one. The characteristic polynomial is J-(X) = X N C(X), where 

N K 

= E ^ L + A- L -2) + 5^j Cj . 

L=l j=l 
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We will find a stable parameterization P corresponding to the fixed point solution 0. Suppose that P 
is of the form P(z) = Ylh=o z ^k- We set, therefore, Po = 0. Also, we chose a solution A of £(A) = 0, 
which amounts to choosing the stable manifold we want to study, and set Pi so that the numerical error 



is minimized. When n > 2, matching coefficients of z n in (35), we obtain 

K 

C(X n+1 )P n+1 + 5 Y, c j sin (i p 



<n\ 



For a generic set of values of 71, 



,7tv and C 1; . 



0. 



(36) 



n+l 



,Ck, we have that £(A n+1 ) / 0, for all n £ N 



Therefore, we can solve the equation (|36|) and get a non-resonant eigenvalue. 

6.2 and assume that we know P— ( n_1 ) and the 



We keep the polynomials P, S J , C J as in Section 



S^Ci corresponding to P-( n x ) up to order n. We use (36) to compute P n and then (31) to compute 



the , C J corresponding to n up to order n — 1. The only difference with the short range case is that, 
when solving the recursion, we need to divide by a slightly different factor. 




Figure 1: Four parameterizations for the Frenkel-Kontorova model of example |6.2| with parameters 
given in table [TJ 



Example 6.2. Consider an specific singular limit. Let the Frenkel-Kontorova model with iV = 3, 
K = 1, 5 = 0.4, and C\ = 1. We fix the values of 71, 72, 73 in four examples. 

Using the algorithm proposed of this section, we find the first 100 coefficients of the Taylor series 
expansion of the parameterizations corresponding to the following values. In each case, the computed 
eigenvalue A corresponds to the slow manifold. 



The solution of the parameterization problem (35) is, in fact, a family of functions that depend on 
the size of the derivative. We find uniqueness only when the derivative P'(0) is fixed. This value is 
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Figure 3: These graphs illustrate the absolute value |3>(.P)| of the error function for the Frenkel- 
Kontorova model of example 6.2 with parameters given in table [I] Logarithmic scale is used in the 
vertical axis. 
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Parameterization 


7i 


72 


73 


5 


P'(0) 


A 


Pa 


1 


0.1 


0.00 


0.4 


10.0 


0.592583231399561 


Pb 


1 


0.14 


0.00 


0.4 


10.0 


0.609158827181520 


Pc 


1 


0.1 


0.01 


0.4 


10.0 


0.603202338024902 


Pd 


1 


0.1 


0.03 


0.4 


10.0 


0.621569001269222 



Table 1: Parameters for the Frenkel-Kontorova examples. 



determined so that the coefficients are of order 1. As input we use the parameters given in table [T] and, 
as output, we get the approximation p^ 100 . 

In the example, parameterizations of the slow manifold are computed. The dimension of the problem 
changes but the method allows the continuation of the solution through the singularity. In other words, 
if we regard the difference equation as a dynamical system, then cases a and b would be maps in M. 
and cases c and d would be maps in M 6 . This collapse in the dimension is problematic if one uses a 
dynamical system point of view, but is manageable when the Lagrangian point of view is used. 

We notice that the four parameterization functions P a , p,, P c , and Pd are similar for small values 
of z. However, the difference equation is singular for the parameters of the first and second examples. 
The numerical results are illustrated in Figure [l| In Figures [2] and [3j we can see an approximation to 
the value of <&(P) near z = 0. In each case, we provide a graph of <I>(P- 100 ). These graphs quantify 
the error in the approximation. 



6.4 The Heisenberg XY model 

Consider the difference equation mentioned in the introduction and given by Q. The characteristic 
polynomial of the fixed point 8* = is P(A) = A 2 — (2 + e)A + 1. The corresponding parameterization 
equations can be written as 

sin(P(Az) - P(z)) + sin(P(A~ 1 z) - P{z)) - e sin P(z) = 0, (37) 

where A is a stable root of T . 



Equating terms of order n in (37) we obtain that, for n = 0, the choice of Pq = corresponds to 
choosing the fixed point solution we are studying. The term of order n = 1 amounts to choosing the 
manifold and setting the numerical scale at which we are working. 

The equations obtained matching order n > 2 are 

P(A" +1 )P n+1 = [sin(P^) -P^(Xz))] n+1 (38) 

+ [sin (P^(z) - P^CA" 1 *))] n+1 + e [sin (P^ "(*))] n+1 • 

The only difference with the Chirikov model or standard map is that we also have to compute S{z) = 
sin(P(z) - P{\z)) and C(z) = cos(P(z) - P(Xz)). Of course, ^A" 1 ^ = sin(P(A- 1 z) - P{z)) and 
C{X~ l z) = cos(P(A -1 z) — P(z)). In the inductive step we assume that we know P^™- 1 ) and the 
S,C,S,C corresponding to P^(™ -1 ) up to order n. Using (38), we can compute P n and then use (31) 
to compute S, C, S, C corresponding to P- n up to order n + 1 and the induction can continue. 
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6.5 Non resonant invariant manifolds in Froeschle maps 

The Froeschle map is a popular model [Fro72| IOV94| of higher dimensional twist maps. It is designed 
to be a model of the behavior of a double resonance. In the Lagrangian formulation, the Lagrangian of 
the model is a function S : R 2 x M. 2 — > R that is given by 

s(e ,e 1 )= l -(o l -e 1 ) 2 + w{e ), 

where W : R 2 — > R is a potential functions such that W(9 + k) = W(9), for all k £ Z 2 . The resulting 
Euler-Lagrange equations are given by 

6 k+1 - 29 k + e k ^ + VW(9 k ) = 0. (39) 

If VW(0) = 0, then 9* = is a fixed point solution and the characteristic function is 

£(A) = det ((A + A~ x - 2)1 + D 2 W(0)) . 

The particular example used by Froeschle is 

W(xx, x 2 ) = acos(27rxi) + &cos(27rx 2 ) + ccos(27r(:ri — x 2 )). 

The matrix I— ^D 2 W(0) is a 2 x 2 symmetric matrix. Typically, it has two real eigenvalues u\ and w 2 . 
It turns out that there are four roots of £(A), that constitute the spectrum of the fixed point. They 
are given by the solutions of 

A + A -1 - 2u)i = 0. 

It is easy to see that these four solutions are given by uii db yjcof — 1, for i = 1,2. From this, we conclude 
that the solutions are of the form Ai,A2, A^ , where |Ai| < | A2 1 < 1. We have to consider three 
possibilities: 

a) < Ai < A 2 < 1. b) < Ai < 1, |A 2 | = 1. c) |Ai| = |A 2 | = 1. 

The classical theory of invariant manifolds allows to associate an invariant one dimensional manifold 
with Ai in cases |aj and|b]), and a two dimensional invariant manifold in case [a]). See also |dlL97| for 
the case < Ai = A 2 < 1. 

The parameterization method allows also to make sense of each manifold tangent to the space 
corresponding to A 2 provided 7^ A 2 and A| 7^ Ax, for k > 2,k G N. The calculations are remarkably 
similar to those of the stable manifolds for the Chirikov map studied before. 

We now indicate the algorithm. Since we are considering the fixed point 9* = 0, we make Po = 0. 
As before, the first coefficient of the parameterization is an eigenvector associated with the eigenvalue 
A 2 . Again, we note that the size of P\ corresponds to different scales of the parameterization, and does 
not affect the mathematical considerations. On the other hand, choosing an appropriate scale is crucial 
in order to minimize round-off errors. 

When W is a trigonometric polynomial, as in the original Froeschle model, we can compute the 
components of \S7W {P~ n ~ 1 )] n ■ We conclude that n— th coefficient of the parameterization satisfies 

((A n + A" n - 2) I + D 2 W(0)) P n = - [VW {P- n ' 1 )] n ■ (40) 
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7 The continuously different iable case 



The approach to invariant manifolds of this paper, also applies when the map Z is finite differentiable. 
Many of the techniques presented above can be adapted to the finite differentiable case. Here, we will 



just present an analogue of Theorem |4.1[ but we leave to the reader other issues such as singular limits 
or dependence on parameters. 

Theorem 7.1. Assume that Z is C r+l . Let A = diag(Ai, . . . , A m ), where A = (Ai, . . . , X m ) is a stable 
non-resonant vector. Assume that Bq is invertible and that r is sufficiently large, depending only on 



A. Then, we can find a C r P solving. (10), P(0) = and such that the range of P\ is the space 
corresponding to the eigenvalues in A. 



The conditions on r assumed in Theorem 7.1 will be made explicit in Lemma 7.3 It will be clear 
from the proof of Theorem |7. 1 1 that , if we assume more regularity in Z, we can obtain differentiability 
with respect to parameters (keeping A fixed). 

We still study the equation &{P) = by implicit function methods but now P ranges over a space of 
finite differentiable functions. To apply the implicit function theorem, we just need the differentiability 
of the functional $ and to show that -D3>(0) is invertible with bounded inverse. One complication is 
that we cannot show that D&(0) is boundedly invertible by matching powers. Nevertheless, we will 
show that -D$(0) is invertible for for functions that vanish at high order. Hence, we will use that we 
can find power series solutions so that we can write P = P < + P— where P < , is a polynomial we will 
assume known and P- are functions vanishing to high order. 

The functional analysis properties of $ will be taken care of in the following lemmas. 

Lemma 7.2. If Z £ (jr+e r ^ ^ ^ then there exists an open subset U C C r {D,W i ) such that the 
operator : IA — > C r (D,M> d ) is C . Furthermore, if£> 1, we have 

N 

D$(P)ip = ^2 djZ (P, P o A, . . . , P o A^) <p o A J '. (41) 



In particular, 



where, as before, B>i = diZ(0, . . . , 0). 



D$(0)<p = Bj (if o A J ) (42) 

j=0 



The Lemma 7.2 is a direct consequence of the results of [dlL099] on the composition operator. Of 
course, the formula (42) is easy to guess heuristically by noticing that it is the leading term in changes 
in P. We emphasize that we are considering A fixed. The differentiability properties of the operator 
with respect to A are much more subtle. 

Let D = W n {p), be the closed disk, in M m , around the origin of radius p. As before, || • ||oo is the 



uniform norm defined in (45). To study C r spaces, it is convenient to recall that higher-order derivatives 
are multilinear maps. Given a C r function g : W 71 — > M d and k < r, then the k— th derivative D k g(a) 
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is a A;— linear map D k g(a) 



symmetric under permutation of the order of the 



arguments. Recall also that there is a natural norm on the vector space of k— linear maps B given by 

||P|| = sup{||fl(fi, . . . ,6)1100 : HCllloc = • • • = Uk\\oc = 1}- 
Once we fix a system of coordinates, D k g can be expressed in terms of the partial derivatives 



1 -D k g(a)(h,...,h)= £ ^ a g{a)h a 



\a\=k 



The usual norm in C r (L',]R ci ) spaces is ||<?||c r = maxj< r sup xgZ) ||P J g(a;)||. 
To study our problem, we define the closed subspace of C r (D, M. d ). 

H r = jp G C r (D,R d ) : D*P(0) = < i < r - l} . 

Clearly P G H r if and only if d a P(0) = 0, for all \a\ < r — 1. The space H r is a Banach space if 
endowed with the norm \\P\\' r = sup {||P r P(:r)|| : x G D}. Because P vanishes with its derivatives at 
and we are considering a bounded domain, this norm is equivalent to the standard C r norm. 



Assume that Bq is non-singular. Then, we can rewrite (42|) as: 

N 



D$(0)tp = B 
where L is the operator defined by: 



V + ^B^Bj&oAi) 



B (Id + L) V , 



(43) 



N 



L(cp) = Y,B 1 B J (tpotf). 



i=l 



Lemma 7.3. Let fi be a real number such that: 
a) ||A|| r <n<l, 



N 



b) ^SWB^BiW <1. 



i=l 



T/ten £/ie linear operator P<I>(0) : P r — >• P r is invertible with bounded inverse. 



We emphasize that if A is a stable non-resonant vector of eigenvalues, the conditions of Lemma 7.3 



are satisfied for sufficiently large r. This is the condition alluded to in Theorem 7.1 



Proof. Since P<I>(0) = B Q 1 (ld + L) it will suffice to prove that L is a contraction in the 
note that if <p G H r , then tp o A J G P r and L</? G H r . In addition, 



| ' r norm. We 



[D r {<p o AP) (x)) (6, ...,&)= [(PV) (A%)] (A^i, . . . , A'£ r )- 



Hence, if x G P then 



|P r (</?o A J ') (x)|| < ||(PV) (A J 'a;)||||A i || r < ||(PV) (A j x)\\ fj, 
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and, since AD C D, we get that \\<p o A^W^ < fjp \\<p\\' r . From this we conclude 

N / N \ 

\\L(<p)\\' r K^Pa'BiUvo P\\' r < [J2^\\ B o lB 4 )M'r- 
j=l \i=l J 

This shows that L is a contraction. □ 



To complete the proof of Theorem 7.1 we argue in a similar manner as in the applications of 



the implicit function theorem before. Following the method in Section 4.5, and fixing Pi to be an 
embedding on the space, we can find a unique polynomial i- >< of degree r — 1 in such a way that 
p<(0) = 0, (P < )'(0) = Pi and Z^'$(P<)(0) = 0, for all < j < r - 1. We recall that the coefficients of 
P < are obtained by matching coefficients. In addition, if we choose Pi small, they will also be small. 
Furthermore, we also note that if P- 6 H r , we have that 

D j $ (P< + P-) (0) = 0; 0<j<r-l. 

The above remarks can be formulated in terms of functional analysis as saying that the operator 
$(P-) = <3?(P < + P-) maps an open subset of H r into H r . 

If we take Pi small without changing its range (this is the same as the change of scale that we 
considered before), we have that $(0) will be as small as desired, D&(0) = D$(P < ) will approach 
D^O) and, in particular will be invertible. The differentiability properties of $ will remain uniformly 
differentiable. Hence, we can deduce Theorem |7.1| from the implicit function theorem in Banach spaces. 

More explicitly, consider the fixed point problem for P- given by: 

-L>$(0)~ 1 $(P < + P-) + P- = P-, 

as acting on H r . When P < is chosen corresponding to a specific small Pi, the left hand side will be a 
contraction and map a ball in H r near the origin into itself. 

8 Final comments 

We have studied invariant objects in the context of difference equations. The goal was to find the analog 
of stable and unstable manifolds that are associated to fixed point solutions of hyperbolic type. We 
generalized the method proposed in [CFdlL03a, CFdlL03b, CFdlL05] and used the implicit function 
theorem in Banach spaces to perform this task. 

The method is not only theoretical but can be implemented numerically as it was shown in example 



6.2. The method is robust in the sense that it can deal with certain singular limits. In particular, slow 
manifolds can be detected and approximated even in the presence of a singularity. 

For future work we plan to consider the finite differentiable case. Also we can use the tools to 
consider the case of conformally symplectic maps and applications to economics. In addition, we 
showed that the variational theory has a smooth dependence on parameters and is robust. We plan to 
explore variational formulations of Melnikov's theory. 
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A Banach spaces of analytic functions 



In this appendix we study spaces of analytic functions (taking values and having range in Banach 
spaces) and, in particular, the composition operator Cf(g) = fog between analytic functions. We show 
that the composition operator is itself analytic when defined in spaces of analytic functions. 

We call attention to the paper |Mcy75| which carried out a similar study and showed that the 
operator T(f,g) = f o g was C°°. The paper |Mey75| showed also that many problems in the theory 
of dynamical systems -invariant manifolds, limit cycles, conjugacies- could be reduced to problems 
involving the composition operator. Using the result of analyticity presented here, most of the regularity 
results in |Mey75| can be improved from C°° to analytic. 

There are of course, specialized books which contain much more material than we need. For example 
|Hof881lNac69] . 

A.l Analytic functions in general Banach spaces 

Let E and F be Banach spaces. We define S k (E, F) as the linear space of bounded symmetric /c— linear 
functions from E to F. For each a k G S k (E,F), the notation a k (x® fc ) denotes the A;— homogeneous 
function E — >• F given by 

afc = a k (x,. . . ,x) . 

On S k (E, F), we require to have a norm || • \\s k (E,F) such that 

\\a k (xi, . . .,x k )\\ F < \\a k \\ Sk (E,F)\\xi\\E ■ ■ ■ \\x k \\E, (44) 
for all a k £ S k (E,F). In particular, this implies that 

a k (x® k ^ < \\ak\\s k (E,F)\\x\\ E - 
Remark A.l. Of course, a natural choice of norms is 

||a fc || = sup{||a fc (xi, . . .,x h )\\ s : \\xi\\ E = ■■■ = \\x k \\ E = 1}, 
but there are others. In the specific case E = C £ and F = C d , we have a found that the norm defined 



in (48) is useful and well adapted to analytic functions of complex variables. 

Once the norms on the spaces S k (E,F) are fixed, we can define analytic functions from E to F. 
Throughout the rest of the appendix, E(5) = {\\x\\ E < 5} will denote the closed ball in E centered at 
the origin with radius 5. 

Definition A.l. Let 5 > 0. The space of analytic functions from E to F with radius of convergence 5 
is the set of functions A$(E,F) given by 



A S (E, F) = If : E(5) -> F f{x) = £ a k (x® k ) , a k E S k (E, F), \\f\\ s < oo , 
I fc=0 J 

where \\f\\s '■= Y^=o \\ a k\\s k (E,F) $ • The a norm \\ • makes A$(E,F) into a Banach space. Given 
p > 0, we will denote A$(E,F) the closed ball centered at the origin with radius p. 
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It is possible to check that if g € A P & (E, F) and / G A p (F, G), then / o g e Ag(E, G). In fact, it is 
proved in |Mey75| that the function T : A P S (E,F) x A p (F,G) -)• given by T(g,f) = f o 9 is 

C°°. In particular, this is true if we consider the Banach spaces of analytic functions that are used in 
the text. Further details are given below for some specific situations that are relevant to our problem. 



In subsection A. 3, we will show that an operator similar to T is analytic. 



A. 2 Analytic functions of complex variables 

We are interested in some concrete Banach spaces of analytic functions. We will define specific norms 
for the following spaces of k— linear functions: 

. S k (C e ,C d ), 

• S k (E, F) where E = A s (C m , C e ) and F = A 5 {C m , C d ). 
First, consider the Euclidean spaces C r with uniform norm 

IMloo = Hfai, • • • >*7r)[|oo = max{|?7i|, . . . , |r/ r |}, (45) 
for each 7/ = (771, . . . , rj r ) E C r . We notice that this norm satisfies 

\z a \ < (IN|oo) W , (46) 

for all multi-indices a 6 W + , and z £ C r . 

Let S k (C , C d ) be the space of bounded symmetric A;— linear functions from C £ to C d . Clearly, the 
spaces Sk(C^ ,C d ) consist of polynomial functions in £ complex variables and d coordinates. It is well 
known that S k (C , C d ) and the space 1-Lk(C e ,C d ) of homogeneous polynomials of degree k are linearly 
isomorphic. 

In particular, given q 6 S k {C l , C d ), the polynomial q (z® fc ) is an element of 'H/ C (C^, C d ). Conversely, 
for each homogeneous polynomial g € %fc(C^,C d ), there exists an unique q € 5fc(C £ ,C a! ) such that 



q (z® fc ) = g{z). Below, in the proof of Lemma A. 2 we will describe how such q can be found 



We need to provide each Sk(C e ,C d ) with a norm that satisfies condition (44). Consider q € 
Sfc(C , C d ). Clearly, the homogeneous function q (z® k ) is of the form 



q(z® k 



£ (47) 

\a\=k 



where ry a are constant vectors of coefficients in C d . Using these coefficients, we define the norm of as 

h\\s k (c e ,c d ) '■= ^2 W^ a W°°- ( 48 

\a\=k 

If it is clear from the context, we will write the previous norm as simply \\q\\. 
Lemma A. 2. The norm on Sk(C^,C d ) defined in (48) satisfies condition (44). 
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Proof. Fixing k > 0, we consider a multi-index a G with \a\ = k. We will give an useful description 
of the unique multi-linear symmetric polynomial function K a G 5fc(C^,C) that satisfies K Q ,(z® fc ) = z a . 

Let 11^ be the permutation group of k elements. Let fa, . . . , fa G be multi- indices such that 
/3i + • • • + fa = a and \fa\ = 1. Using these vectors, we define the function n a G 5fc(C , C) as 



Ka^l, . . . , Z k ) :— ^ Z n(l) Z t(2) ' ' ' Z w(k)> 



(49) 



where Zi,...,z k G C. It is clear that K a is a well-defined multi-linear function and K a (z® fc ) 
^(z, 2, . . . , z) = z a . In addition, using (46), we get that K a satisfies 

\K a (zi, ... ,Z k )\ <— ^ 1^(1)11^(2)1 "' '" 



7ren fc 



1 „ 

-k\2^ ^ 



) ||oo |p7r(2) lloo " ' ' 1 1 Z-ir(k) I 1 oo — \\ z l\ 



\ z k\ 



7ren fc 



From this, we conclude that if q G <Sfc(C , C ) and q satisfies (47) then 

q (zi, . . . , z k ) = K a (zi, ZkjVa, 

\a\=k 



and therefore 



j Zk) | |oo 5; ^ ] ■ • • j -2fe)| || ^?a||oo — l|(?lls fc (C f ,C d ) 1 1 z l I loo " " " ||^fc||oo- 

\a\=k 



The last inequality implies that the norm defined in (48) satisfies (44). 

From the discussion above, we also conclude that if g G A P (C £ , C d ) then g is of the form 

oo 

g(z) = zClr]a - 

k=0 \a\=k 



□ 



(50) 



Therefore, with the norms defined in (48) on the space of k— symmetric functions Sk(C^,C d ), we get 



the following norm, for any g G AJC , C 



i r d\. 



If Hp = £ 2 H^ll 00 

fc=0 \|a|=fc 



P 



(51) 



This norm has the following property. If g G A p (C e , C) and r\ G C rf then g • 77 G j4 p (C^, C rf ) and 

ll^-^llp = llffllplhlloo- (52) 

We also have the following. 

Lemma A. 3. If g\,g2 G A p (C r ,C), then g\g2 G yl p (C r ,C) and 

Hfllfl^Hp < HfflllplMlp- 
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Proof. Let 51,52 £ ^4p(C r ,C). Then, they are of the form 



«(*) = E E ^a. 
fc=0 \a\=k 



for i = 1,2. This implies that 



9l (z)g 2 (z) = E E E za+ ^4 

ii,t2=o\ct\=h \P\=t* 

00 00 / \ 

= E E ^ = EE'1 E « > 

^=0 |q| + |/3|=£ ^=0 | 7 |=f \a+^=7 / 

for all z £ C r (p). Notice that the coefficients rf a are complex numbers and therefore they satisfy 



Wo 



. Using (51), we get that 



isis2ii P = E E 

£=0 \y\=i 



El 2 
VaV/3 

o+/3=7 



/<EE E 

£=0 | 7 |=£a+/3=7 



E E kV 1 EE \4\p £ ' l wimutm,,- 

^i=0|o|=4 / \^2=0 |/3|=£i 

This also shows that the series involved in || 31 32 Hp converges. 

The previous result also shows that A p (C r ,C) is a Banach algebra. Furthermore, if g £ A p (C r ,C s 
and a £ is a multi-index, then g a £ A p (C r , C) and 



□ 



| «|| < || fl ||H 
\y lip ^ \\y\\p ■ 



(53) 



A. 3 The composition operator 

Let E = A 5 {C m ,C £ ) and F = A 5 {C m ,C d ) be two Banach spaces of analytic functions. We will 
define norms on the corresponding spaces of A;— symmetric linear and bounded functions. For each 
bk £ Sk (E, F) we define the norm 



\\h\\s k (E,F) '•=sup{\\b k (g 1 ,...,g k )\\s : \\gi\\s 



\9k\\s = !}• 



(54) 



Notice that this norm satisfies condition (44). With these norms, we can show that each analytic 
function in A P (C £ , C d ) can be extended to an analytic operator in A P (E, F), acting on spaces of analytic 
functions. 

Lemma A.4. Let f £ A p (C E , C d ) and 5 > 0. Let E = A s (C m , C £ ) and F = A 5 (C m , C d ) be two Banach 
spaces of analytic functions. If Cf : E(p) —> F is the operator given by Cj(g) = f g, then Cj is well 
defined and analytic^ In addition, with the norms on the spaces Sk(E,F) defined in (54), we have that 

\\Cf\U p (E,F) < ll/llp- 



J This means that Cf £ A P (E,F). 
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Proof. We know that / is of the form 



f( z ) = ^2 a k 



k=0 



with a k G S k {C\C d ), and 

oo 

II/IIp = ^2 hk\\s h {&,&) P k < °°- 
k=0 

Each multi- linear function a k G tSfc(C^,C rf ) can be extended to a function a& G 5fc (E,F) by 

[ajfeCffi, • • .,5fc)] 00 : = a k (gi(z), . . .,gk(z)), 
for each z G C m such that ||z||oo < and all gi, . . . , g k G -B. 



As in the proof of Lemma A. 2 we can assume that each a k is of the form 



a k (zi, . . . , z k ) = K a (zi, • • • ,^)r? a - 

|a|=fc 

In particular, each polynomial K a can be extended to a function k a G 5^ (E,F). This implies 

afc ■■■,9k) = ^2 Ra (9i, ■ ■ -,gk)i]a- 

\a\=k 



In addition, using (49) and inequality (53), we get that 



Using (52) 



\& a (gi, ■ ■ -,gk)\\s < \\gi\\s ■ ■ ■ \\gkh- 



\a k (gi, ■ ■ ■ ,gk) h < ( ^ l r /a Joe I !<> ■ ■ ■ 1 1. <y/.- 1 U- 

Aa\=k 



Therefore, for all k > 0, the norm (54) on S k (E,F) is such that 



ll«fclU fc (£;,F) < \\ a k\\s k (C e ,C d )- 
Finally, we notice that the operator Cf can be written as 



fc=0 |a[=Jfc 



and conclude that 



\Cf\\ Ap (E,F) = Yl \\*k\\s h (E,F) P k <Y1 W a k\\s k 



(C e ,C d ) P 



< 00. 



k=0 



k=0 
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